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Abstract 


This paper proposes a numerical approach to solve a two-layer shal- 
low water formula with a slope and uneven bottom. The finite volume 
method (FVM) is applied to solve the shallow water model because the 
method is suitable for computational fluid dynamics problems. Rather 
than pointwise approximations at grid points, the FVM breaks the domain 
into grid cells and approximates the total integral over grid cells. The 
shallow water model is examined in two cases, the shallow water model in 
the steady state and the unsteady state. The quadratic upstream interpo- 
lation for convective kinetics (QUICK) is chosen to get the discretization 
of the space domain since it is a third-order scheme, which provides good 
accuracy, and this scheme proves its numerical stability. The advantage of 
the QUICK method is that the main coefficients are positive and satisfy 
the requirements for conservativeness, boundedness, and transportation. 
An explicit scheme is used to get the discretization of the time domain. 
Finally, the numerical solution of the steady state model shows that the 
flow remains unchanged. An unsteady-state numerical solution produces 
instability (wavy at the bottom layer). Moreover, the larger slope results 
in higher velocity and higher depth at the second layer. 
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1 Introduction 


A flow with constant density where the horizontal scale is much bigger than 
the vertical scale is called shallow water. Shallow water is a thin layer of 
fluid of constant density in hydrostatic equilibrium, surrounded by a solid 
surface below and a free surface above, and is likely in contact with another 
low-inertia fluid. The simplest and most useful models in geophysical fluid 
dynamics are single-layer models because they can explain the rotational ef- 
fect with a simple structure without considering the effects of complex strat- 
ification. We can investigate the effect of the layers by adding a layer to 
them. The two-layer in the geometry of shallow water is not only a layered 
liquid model, but also it is an excellent model of a wide range of ocean and 
atmospheric phenomena; see [23]. 

One way to depict fluid flow is as a multilayer flow, where one layer 
flows on top of another. This difference in a layer is caused by temperature 
differences, which lead to different densities. Ocean water can be viewed as a 
fluid with two layers. The temperature of the water in the first layer, which is 
sunlight-exposed water, is greater than the temperature of the second layer, 
which is below the first layer; see [23]. 

Shallow water flow in two layers, especially its exact and numerical so- 
lution, is an interesting subject to be studied. They particularly find the 
best-governing equation that can describe physical phenomena. 

The two-layer fluid flow model over the inclined and flat bottom was 
investigated in [17]. The model utilized a non-Newtonian fluid. They used 
the flow of mud in the water to illustrate the model. Then, this model was 
solved analytically and numerically. The McCormack’s scheme was used to 
obtain the numerical solution since this scheme is suitable for systems that 
meet the hyperbolic conservation law, and it has second-order accuracy. The 
numerical simulation of stratified fluid, supposed to be two-layer shallow 
water flows through a channel with irregular geometry, was studied in [5]. 
They proposed extensions of the Q-schemes of van Leer and Roe to get the 
numerical solution of the model. The central upwind was used to solve the 
two-layer shallow water equations [10]. 

Furthermore, the fluid flow model with inclined and uneven bottom was 
developed in [25]. In this model, it was assumed that the fluid had a single 
layer and flowed in horizontal and vertical directions. The constructed model 
was completed analytically and numerically using the finite difference method 
by applying periodic boundary conditions. A numerical scheme stability 
analysis was also performed using the Fourier analysis developed by Von 
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Neumann. The stability analysis revealed that the flow was stable, which 
means that the error of a numerical solution will not grow infinitely with 
time. 


The analytical and numerical methods were developed quickly to find the 
best solution for the model. The analytical method is the best method to find 
the solution of the mathematical model since this method does not have any 
errors that appear from the truncation of the method itself. Unfortunately, 
analytical methods sometimes require some complex calculations and are very 
difficult to solve. To handle this problem, numerical methods can be used 
to approximate the solutions of the governing equations. The classical finite 
difference method, in which derivatives are approximated by finite difference, 
is a well-known method to approximate the solution. Discontinuities become 
a problem when they occur in the solution. Such a method can be used as 
well with certain treatments near discontinuities in the solution. The finite 
volume method (FVM) is a method based on an integral form instead of the 
differential equation. Rather than pointwise approximations at grid points, 
the FVM breaks the domain into grid cells and approximates the total integral 
over grid cells. This integral is divided by the volume of the cell. These values 
are changed at each step by the flux passing through the edges of the grid 
cells. The main problem is to determine the proper numerical flux function 
that properly approximates the correct flux [11]. 

The study of the two-layer shallow water system (flat bottom) analyt- 
ically and numerically was investigated in [3]. For an analytical solution, 
they considered a relaxation approach, which offers greater decoupling and 
accessible eigenstructures. The numerical solution of a two-layer shallow wa- 
ter equation using the FVM with an explicit scheme was investigated in [7]. 
Then, Swartenbroekx, Zech, and Soares-Frazao [22] predicted the bed load 
transport induced by dam break waves by proposing two-dimensional (2D) 
two-layer shallow water equations, where the upper layer was made of clear 
water, and the lower layer was made of a dense mixture of water and moving 
grains. The governing equation of the model was solved by applying Harten— 
Lax—Van Lear Finite Volume Scheme and a well-balanced Roe method [9]. 
Recently, Muhammad [13] simulated the flowing fluid using the FVM via 
OpenFOAM. 

On the other hand, the flat-bottomed shallow hydraulic analogy of the 
Gibraltar Strait was developed with wind stress in mind [6]. In this model, 
the flow was assumed to be quasi-horizontal. The y component was ignored 
,and the density was uniform in each layer. The model was numerically 
solved using the FVM based on the Roe-Riemann scheme approach. Chi- 
apolino and dan Saure [8] also derived a two-layer shallow water flow model 
with some limitations by ignoring the vertical velocity component and as- 
suming that the velocity is uniform across the cross-sections of each layer. 
This model was solved numerically using the Godunov scheme, which is an 
approximation of the Riemann problem at each cell boundary. A method 
for the numerical solution of the nonlinear regularized long wave equation 
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was studied in [20]. They [19] proposed an implicit method for the time 
discretization based on the 6-weighted and finite difference methods, while 
the spatial discretization is described with the help of the finite difference 
scheme derived from the local radial basis function method (the local collo- 
cation method). The localized meshless collocation method was studied in 
(16, 15]. A new weighted essentially non-oscillatory (WENO), a finite dif- 
ference Hermite RBF-WENO scheme, and a modified weighted essentially 
non-oscillatory (WENO) scheme are procedures for solving hyperbolic con- 
servation laws were studied in [18, 1, 2]. 


The study conducted in [4] obtained a numerical solution of the advection- 
reaction-diffusion equation using the FVM. In the FVM, the integral form 
is used to discretize the equation. The discretization of space used the 
Quadratic Upwind Interpolation for Convective Kinematics (QUICK) scheme, 
while the discretization of time used an explicit, implicit, or Crank Nicolson 
scheme. 


In this paper, we solve numerically a mathematical formula of a two- 
layer shallow water flow with slope and uneven bottom by using the FVM 
with a QUICK scheme to discretize the space domain and explicit scheme 
to discretize the time domain. The new one in this article is how to find a 
solution of the model when the bottom of the shallow water is not flat. The 
FVM is used to solve the model because this method is suitable for compu- 
tational fluid dynamics problems. The QUICK method is chosen since it is a 
third-order scheme that provides good accuracy, and this scheme have been 
proven its numerical stability. The advantages of this method are that the 
main coefficients are positive and satisfy the requirements for conservative- 
ness, boundedness, and transportiveness [24]. A numerical simulation of the 
solutions is presented to describe the behavior of this model. 


The organization of this work is as follows: Section 2 derives the governing 
equations of the models in steady and unsteady states. In Section 3, the 
FVM is constructed as a numerical solution with the QUICK scheme for 
discretization. Section 4 provides numerical simulations of the shallow water 
flows with several given parameters. Finally, Section 5 ends up with a brief 
conclusion. 


2 Governing equations 


The governing equation of the two-layer shallow water with slope and the 
uneven bottom is described in Figure 1. It is a non-dimensional form after 
deriving from continuity and momentum equations of shallow water equa- 
tions. For detailed calculating, revisit [12]. The equation reads as follows: 


° (m1 — 2) + ° (ui(m — n2)) = 0, (1) 


ax 
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On2 ) 


b)) =0 2 
at + Dy (u2lna — 8) = 0. (2) 

0 0 (a) 

7 + uy a = —gcos6 ae +gsin6, (3) 

0 0 0 (6) 

FP + U2 a = (g' g) cos 0 — g! cos 0 + gsin0@. (4) 
We introduce the non-dimensional variables as follows: 

- D~ 2 
x= La, z= Hz, ome Fae m = Am, 


m= Hip, b=Hb, wm =UHG, @=UHe, 
and 


1 U2 g 
U = (Hg’)2, L=—, YY S95 
g g 
where U is a typical velocity, L is typical length, and T is a typical time scale, 
and H is mean height. The two-layer shallow water with slope and uneven 
bottom in non-dimensional form is 


Om One Ourm OuyN2 


at Ot |SCOOxt an 2) 
Ong Ou2N2 Ougb -_ (6) 
Ot Ox Ox 

du. O (ut\  cosOOm 

at Oa ( 2 ) «Ox Be (7) 
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dur O (uz\ | 1 jm Aetna = 
e+e(2 =(1 : cos 0 cos 0 + sind. (8) 


We call the bottom layer as layer one, and the upper layer as layer two. 
We introduce the variables used in each layer, pi (x,t) and p2(z, t) are pressure 
forces in layers one and two. The upper layer has a pressure force pi(2,t), 
which is lower than the pressure force at the bottom layer po2(z,t). Then, 
m(a,t) is the free-surface movement vertically while 72(az,t) is the distance 
between layers one and two’s interface height, where 72(#,t) < z < 1\(2,t) 
and b(a) < z < m2(a,t). The thickness of the water columns is h,(z,t) 
and ho(x,t). The uneven bottom topography is b(a#), where, in general, 
ho(ax, t) = no(a,t) — b(x), for the flat bottom no(x,t) = he(x,t), and @ is the 
slope of the bottom surface. Moreover, y = g/g, where g is the gravitational 
force, and g’ = g(f i 1) is the reduced gravity, while py and p; are density in 
layers one and two. The proposed model, (1)—(4), is a powerful model of many 
geophysically interesting phenomena as well as being physically realizable in 
the laboratory [23]. 


2.1 Steady-state two-layer shallow water flow 


In this part, shallow water is assumed to be in a steady state. This means at 
every point, the mass fluxes over the depth and energy are constant [14]. It 
is assumed that the height and water debit are constant with respect to time 
[21]. Equations (5)—(8) become 


Ourm _ Ourne 


Ox ax . (9) 
rae Sg wo 
g = en +sin 0, (11) 
0= (1- ~) cos 0 — cos 9 Me +sin 6, (12) 


where x € [0, K] in which K is a constant at the boundary space domain. 
The equation system (9)—(12) is a function of spatial variable, such that the 
initial conditions for the equation system (9)—(12) are given as follows: 


m(a,0) = tan dx + Asin(3z), (13) 
No(x,0) = tan 0x + Asin(3x) — 3, (14) 
b(a) = tan 0x + K tan@ + A(cos(5a) — sin(3a)) — 4, (15) 


and the boundary conditions are given by 
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m (0) = tan @(0) + Asin(3(0)), (16) 
n2(0) = tan 6(0) + Asin(3(0)) — 3, (17) 
m(K) = tan d(x) + Asin(3(K)), (18) 
no(K) = tand(i) + Asin(3(K)) — 3, (19) 

on = ore = tan 67, (20) 


where A is wave amplitude. 

The steady state of (9)-(12) can be applied to river at rest conditions. 
It seems that the river consists of two-layers. The temperature of the first 
layer will be lower than the temperature of the second layer when the river 
is exposed to sunshine. As a result, the first layer pressure is lower than 
the second layer pressure. The layers compete with one another when the 
interface between the two layers varies from position to position (wavy). 
Similarly, on an uneven base, generally, topography and bottom layer will 
exert forces on one another. Because it is in a steady state, the bottom 
shape is unevenly followed by the first and second levels. We will show this 
in numerical simulations. 


2.2 Unsteady-state two-layer shallow water flow 


For an unsteady state two-layer shallow water flow, it is assumed that the 
flow of shallow water at the upper free-surface, the gravitational force, is 
imposed on the fluid as the restoring force. Because we consider hydrostatic 
balance, the changes of surface velocity, and height are quite small around 
the depth of the h surface wave at constant conditions [23]. The depth of 
the fluid is as much as the displacement or movement of the free surface. We 
linearize (5)-(8) by defining new variables to describe the displacement of 
the fluid at rest; that is, 


hy(a,t) = Ay + hy(a,t) = Ay + (ni (2, t) — m9(z, t)), (21) 
ho(ax,t) = He + hg(a,t) = He + (ng(a,t) — b(a)), (22) 
ur(z,t) = uj (2, t), (23) 
(0) — tat (24) 


Equations (21)—(24) are substituted into (5)-(8). Then, the variable fol- 
lowing the derivation of steady state two-layer shallow water flow is non- 
dimensionalized. We obtain the following equations: 


On One _ 

a a” 25) 
Ong Oug Ougb _ 

i Op Ge ee 
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Ou. cos0 Om 


ind 2 
Ot y O ae 7) 
Ou 1 On, One ; 

je omy (ts 56 —¢o8 O24 2 
Fi (1 - | cos 6 Ae 0089S, + sin 0, (28) 


where x € [0, K],¢ € [0, L], with the initial conditions of the flow is stationary 
and the speed of the water remains along the channel 


m(a,0) = tandx + Asin(3z), (29) 
no(a,0) = tan dx + Asin(3x) — 3, (30) 
b(x) = tan Ox + K tan 6 + A(cos(5x) — sin(3x)) — 4, (31) 
u1(x,0) = ue(x, 0) = 0, (32) 


and the boundary condition is assumed that the flow is not moving, so that 
the velocity is equal to zero. The boundary condition of unsteady state two- 
layer shallow water flow is periodic functions as follows: 


m(0,t) = tan 6(0) + Asin(3(0)), (33) 
no(0,t) = tan 6(0) + Asin(3(0)) — 3, (34) 
m(K,t) = tand(K) + Asin(3(4)), (35) 
no(K, t) = tan 0(K) + Asin(3(K)) — 3, (36) 
om (O.8) = onlO.4) = tan 67 + 3Acos(3(0)), (37) 
uz (0,t) = ua(0,t) =0, (38) 


where x > 0 andt > 0. 


3 Finite volume method (FVM) 


The FVM is a technique using integral form rather than differential equations. 
Pointwise estimates at grid points are preferred. The domain is divided into 
grid cells by the FVM, which then approximates the entire integral over the 
grid cells. This integral is divided by the volume of the cells. The value 
changes at each step depending on the flow passing through the edges of the 
grid cells. The main problem is determining the appropriate numerical flow 
function that is directly close to the true flow [11]. 


By following the procedure in [24], we integrate the governing equations 
for steady and unsteady state, (9)-(12) and (25)—(28), respectively, along 
control volume and time. Furthermore, we use QUICK scheme to discretize 
space domain and explicit scheme to discretize time domain. 


Here, we summarize the procedure of the proposed method as follows: 
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1. Grid generation 

The first step in the FVM is to divide the domain into discrete control 
volumes. Let us place a number of nodal points in the space between 
A and B. The boundaries (or faces) of control volumes are positioned 
mid-way between adjacent nodes. Thus each node is surrounded by a 
control volume or cell. It is common practice to set up control volumes 
near the edge of the domain in such a way that the physical boundaries 
coincide with the control volume boundaries. 


2. Discretisation 
The key step of the FVM is the integration of the governing equation 
(or equations) over a control volume to yield a discretized equation at 
its nodal point P. 


3. Solution of equations 
The resulting system of linear algebraic equations is then solved to 
obtain the distribution of the solution at nodal points. 


3.1 Steady-state FVM 


To solve (9)—(12) by using the FVM, first, we integrate (10) with respect to 
Control Volume (CV) as follows: 


Ou2N2 Ougb = 
pe ao st) av =0. (39) 


This is the equation related to the first layer of shallow water flow. Since the 
model of two-layer shallow water flow has one dimension in space, the CV 
can be described as in [24]. Then (39) becomes 


A (U2 _ Uw) (N26 Tw) =A (Ure = U2w) (be ~ bw) ’ (40) 


where index e is the east face and w is the west face. By applying QUICK 
discretization, we transform the control face to be nodal and after simplifying 
(40). Then we get 

1 7 3 3 1 7 3 


3 
ginww — ginw + amp 4 gine = gbww gow goP g0E- (41) 


Equation (41) can be written as a linear equation of the system by substi- 
tuting i = 2,3,...,N —1 and applying initial and boundary condition such 
that we have 


Iran. j. numer. anal. optim., Vol. 13, No. 2, 2023,pp 317-335 


Habibah, Lina and Kusumawinahyu 326 


3.63 
ry 3 00---0 0 0 0 
—7 3 
ces 0 0 N22 - 
23 3 
—7 3 
Se ee 24 Fy 
a f Nel en | (42) 
-0 0 0 0 
1-7 3 8 12K FrR-2 
8 8 8, 8 N2K-1 K-1 
0 0 00---0 = —-= 
- § 8 8A 
with 
ie 1 1 ed 
an + ~AgxT1 4 a { Spend - iieh Az tan(6)y 
R ic a ae en 
2 bj bj-1 + -—b; + -b; = 
P; | 8 4 a 8 4 : 8 4 38 i+1 3 2% 2 
4 qm geici ge goett 
FrR-2 1 7 3 3 
Fics bi—2 bi-1 + bi + Sbi41 
1 2 7 2 3 83 e 3 
L gli? gent + ge + goitt g Mtl | 


where T, = tan 0+ A(5sin(—5(0)) —3 cos(—3z(0))). The methods for solving 
a system of linear equations can be used to determine the value of 1;, for 
example, Thomas’s algorithm. 


Furthermore, similar to the first layer, the discretization equation of the 
second layer is obtained by integrating (9) with respect to CV 


Our, OurN2 _ 
[ (2m 28) avo a 


Discretization QUICK is applied to (43) by following the procedure in [24] 
and applying initial and boundary conditions. Then, we arrive at a system 
of linear equations as follows: 


[ 2 Gino & 8H) 
SF Bin ied 0 0 M12 G2 
7 ms | | & 
—7 3 
eee eee a a ae m4 en 
8 8 88 : _ . ; (44) 
oe ed ie We : 
dk, “ie Ss 1-7 3 3 MK-2 Gr_2 
ep yap lm Lew 
—7 
0 0 00-:---0 = —-= 
a 8 8 8-7 


with 


Iran. j. numer. anal. optim., Vol. 13, No. 2, 2023,pp 317-335 


327 Two-layer shallow water formula with slope and ... 


6 1 3 3 6 1 4 
~ glzi-1 + >Axtan(0)y +4 gi g itl gi2i—1 3 Se tan(0)y 
Cs p 7 ae: 1 
G3 g 2i-2 g2i-1 t 3 7% t g itl g Mi-2 
1 7 3 3 
Ga gl2i-2 — gMi-1 + omit oMit1 
Gr_-2 1 ve 3 3 
Gr_i g 2i-2 gi-1 g l24 g 12i+1 
3 3 
L g 21-2 g 2i-1 t g Mi t g itl g Titi | 


The value of 7;; can be obtained by the methods for solving system of linear 
equations. 


3.2 Unsteady-state FVM 


For unsteady-state two-layer shallow water flow, we integrated (25)—(28) with 
respect to CV and time. The procedure to get the integration (25)—(28) with 
respect to CV in the unsteady state, shallow water flow is similar to the 
procedure to calculate the integral of the governing equation in the steady 
state. We apply QUICK discretization. The difference is we integrate the 
governing equation with respect to time, and apply an explicit scheme for 
discretization. First, we integrate (25) to yield 


t+At 
/ ) (s 2 a dV dt = 0. (45) 
The result of the integration of (45) is 
(mp —mp)AAz — (np — nap) AAx = 0, (46) 


where 7,p and 72,p are the height of surface at node P and time t+ At, ny, P 

and 79, are the height of surface at node P and time t. By simplifying (46) 

and writing the index P becomes i, we get the height of the first layer 
met =m? + met? — mi, (47) 


where the subscribe n denotes time t, and n + 1 denotes time t + At. We 
apply the procedure in [24] for the rest of (26)—(28), so that we obtain 


Ar 8 8 8 

1 7 3 3 
1 Pei eg ee 4 
( € i-2— git. + 8 )) (48) 
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At cos@ (1 7 3 3 
tier =u? (Gnts ga stig + moa) 


An ¥ 8 8 8 
+sin OAt. (49) 
ur" oo Uoe + Ae (1 _ +) cos 6 (Gms g Mi-1 gs 7 zn) 
At 1, — Die ye ees Si ; 
er cos 8 (Jmts git + gli + zi + sin OAt. (50) 


The explicit solver is easily programmed, yet the time step t must be 
limited to stabilize the scheme. Experience has shown that we can simply 
enforce the Von Neumann criteria for the single step algorithm, 


At Z 
Ar 


1 


7 (51) 


4 Numerical simulation 


Numerical simulation aims to illustrate movement of shallow water flows 
based on the numerical scheme obtained in the previous section. Numerical 
simulations are performed under the environment of MATLAB R2018a based 
on MacOS MOJAVE version 10.14.6, which the processor is 1.6 GHz Intel 
Core i5. Numerical solutions of steady-state of two-layer shallow water flow 
are calculated at x € [0,5] by taking Az = 0.1 and y = 1 [17]. Figure 2 
depicts two shallow water layers at rest, where the water velocity is zero and 
there is no time-dependent flow movement. Changes at the bottom surface 
are followed by depth changes at the first and second chamber levels. When 
we increase the slope, the depth of the second layer is higher than the first 
layer. The two-layer shallow water flow conditions remain stable for each 
slope of (0). Negative sign of slope has physical interpretation, that is, water 
flows from the left side to the right side. Figure 3 shows the velocity profiles 
of steady-state of two-layer shallow water flow at t = 0.025. 


Table 1: Comparison of the CPU run times of the proposed method depends on the 
slopes. 


é?=—10 ? = —30 6 = —60 
1,376413 second | 1,536349 second | 1,544937 second 


Furthermore, the numerical simulation for unsteady-state two-layer shal- 
low water flow can be seen in Figures 4 and 5. The parameter (7), slope (9), 
and space domain for the simulation are chosen with the same values for the 
steady-state. The time domain is t € [0,0.1] with At = 0.001. The simulated 
water wave height and velocity are calculated when ¢ = 0.025 and t = 0.1. 
The height profile of unsteady-state two-layer shallow water flow at t = 0.025 
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0=-10° 0 = -60° 


—h1(x,t)) 
—n2(x,t) 
—b(x) 


5 
—n1(x,t) 
—n2(x.t) 0 
— (x) 
5 


c=) S 
D2 2-10 
2 2 
15 
4 
a oe . 
4 1 . x 25 
0 1 2 3 4 5 0 1 2 3 4 5 


Figure 2: The numerical solution of steady-state of two-layer shallow water flow with 
6= 10°, 6= 30°, and @ = 60°. 


0 = -10°, t=0.025 0 = -30°, t=0.025 
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-0.04 
0 


0 = -60°, t=0.025 


v2 


velocity 
2 


0.1 
0 


Figure 3: The velocity profiles of steady-state of two-layer shallow water flow at t = 
0.025. 


is shown in Figure 4, while the velocity profiles of unsteady-state of two-layer 
shallow water flow at t = 0.025 is shown in Figure 5. 


0 = -10°, t=0.025 0 = -30°, t=0.025 


height 


Figure 4: The height profiles of unsteady-state two-layer shallow water flow at t = 0.025. 
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Figure 5: The velocity profiles of unsteady-state of two-layer shallow water flow at 
t = 0.025. 


The numerical solution for unsteady shallow water produces instability 
(see Figure 5). It is because there is a compressive force between the layers. 
Similarly, if the liquid’s bottom is not flat, the bottom layer and topography 
tend to push against one another. This type of force, known as form drag. 
Moreover, in Figure 5, when the slope is increased, the depth of the second 
layer is higher than the first layer, and the second layer has a more wavy 
shape than the first layer. In Figure 5, the flow has bigger velocity when 
the slope of bottom layer is increased. From Table 2, the velocity at layer 
one increases when the slope is increased. Similarly, the velocity at layer 
two increases when the slope is increased see Table 3. The last row in each 
Tables 2 and 3 shows the rate of velocity u,; and u2. The CPU run times of 
the proposed method is written in Table 1. 
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Table 2: Comparison of velocity u; depends on the slopes. 
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u, at9?=-—10 | uz at O=—30 | uw at O = —60 
0,01 0,01 0,01 
0,01 0,01 0,01 

0,010000047 0,010000156 0,010000486 
0,010000042 0,01000014 0,010000435 
0,010000035 0,010000115 0,010000356 
0,010000025 0,010000084 0,010000257 
0,010000015 0,010000049 0,010000148 
0,010000004 0,010000013 0,01000004 
0,009999994 0,00999998 0,009999941 
0,009999986 0,009999953 0,009999861 
0,00999998 0,009999934 0,009999804 
0,009999976 0,009999924 0,009999775 
0,009999977 0,009999924 0,009999776 
0,00999998 0,009999934 0,009999806 
0,009999986 0,009999954 0,009999864 
0,009999995 0,009999982 0,009999946 
0,010000005 0,010000015 0,010000046 
0,010000015 0,010000051 0,010000154 
0,010000026 0,010000086 0,010000263 
0,010000035 0,010000117 0,010000361 
0,010000043 0,010000141 0,010000439 
0,010000047 0,010000157 0,010000489 
0,010000049 0,010000162 0,010000504 
0,010000047 0,010000156 0,010000485 
0,010000042 0,010000139 0,010000431 
0,010000034 0,010000114 0,010000351 
0,010000025 0,010000082 0,010000251 
0,010000014 0,010000047 0,010000142 
0,010000003 0,010000011 0,010000034 
0,009999993 0,009999979 0,009999936 
0,009999985 0,009999952 0,009999857 
0,009999979 0,009999933 0,009999801 
0,009999976 0,009999923 0,009999774 
0,009999977 0,009999924 0,009999776 
0,00999998 0,009999935 0,009999809 
0,009999986 0,009999956 0,009999868 
0,009999995 0,009999984 0,009999951 
0,010000005 0,010000017 0,010000051 
0,010000016 0,010000053 0,01000016 
0,010000027 0,010000088 0,010000268 
0,010000036 0,010000118 0,010000366 
0,010000043 0,010000142 0,010000442 
0,010000047 0,010000157 0,01000049 
0,010000049 0,010000162 0,010000504 
0,010000047 0,010000155 0,010000483 
0,010000042 0,010000138 0,010000428 
0,010000034 0,010000112 0,010000346 
0,010000024 0,01000008 0,010000245 
0,010000013 0,010000045 0,010000136 
0,010000003 0,010000009 0,010000028 
0,01 0,01 0,01 
0,510000685 0,51000228 0,510007164 
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Table 3: Comparison of velocity ug depends on the slopes. 


ug at9=—10 | ug at O= —30 | ug at Od = —60 
0,01 0,01 0,01 
0,009999236 0,009999431 0,009999728 
0,009999233 0,009999349 0,009999664 
0,009999146 0,00999918 0,009999542 
0,009999225 0,009999147 0,009999474 
0,009999545 0,009999439 0,009999622 
0,010000052 0,010000062 0,01000004 
0,010000542 0,01000064 0,010000417 
0,010000821 0,010000851 0,010000502 
0,010000857 0,010000764 0,010000407 
0,010000736 0,010000576 0,010000283 
0,01000055 0,01000039 0,010000182 
0,010000354 0,010000236 0,010000106 
0,010000169 0,010000109 0,010000049 

0,01 0,010000001 0,010000001 
0,009999849 0,009999904 0,009999958 
0,009999721 0,009999815 0,009999917 
0,009999628 0,00999974 0,00999988 
0,009999584 0,009999691 0,009999852 
0,009999598 0,009999681 0,009999842 
0,009999659 0,009999714 0,009999853 
0,009999739 0,00999977 0,009999879 
0,009999804 0,009999822 0,009999904 
0,009999833 0,009999847 0,009999917 
0,009999824 0,009999841 0,009999915 
0,009999795 0,009999816 0,009999902 
0,00999977 0,009999793 0,009999889 
0,009999785 0,009999802 0,009999893 
0,00999987 0,00999988 0,009999935 
0,01000004 0,010000044 0,010000027 
0,010000262 0,010000256 0,010000145 
0,010000465 0,010000426 0,01000023 
0,01000058 0,010000443 0,010000214 
0,010000494 0,010000344 0,010000159 
0,010000353 0,010000229 0,010000102 
0,010000189 0,010000116 0,01000005 
0,010000014 0,010000006 0,010000002 
0,009999831 0,009999893 0,009999953 
0,009999641 0,009999768 0,009999897 
0,009999446 0,00999962 0,009999826 
0,009999266 0,00999945 0,009999736 
0,009999147 0,009999283 0,009999632 
0,009999163 0,009999199 0,009999553 
0,009999371 0,009999324 0,009999592 
0,009999739 0,009999699 0,009999808 
0,010000128 0,010000137 0,010000082 
0,010000397 0,010000409 0,010000239 

0,0100005 0,010000471 0,010000259 
0,010000565 0,010000492 0,010000257 
0,01 0,01 0,01 


0,509997101 


0,509997386 


0,509998561 
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5 Conclusion 


In this research, we offered a numerical solution for a two-layer shallow water 
flow with slope and uneven bottom for both steady-state and unsteady-state. 
The QUICK and explicit schemes were used to determine the numerical so- 
lution through the FVM. The steady-state numerical solution yielded stable 
results, while unsteady-state numerical solutions produced unstable results. 
The momentum that can be added to or taken from a flow, is what leads 
to this instability. This happened when there are differences in the position 
or surface wave shape between the first and second layers. The layers exert 
a compressive force on one another. Similar to this, if the fluid’s bottom is 
not flat, the bottom layer and the topography will generally exert forces on 
one another. Moreover, the simulations yielded the larger slope resulted in 
higher velocity and higher depth at second layer. 
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